module ResistanceCanopyStomataBallBerryMod

!!! Compute canopy stomatal resistance and foliage photosynthesis based on Ball-Berry scheme

  use Machine
  use NoahmpVarType
  use ConstantDefineMod

  implicit none

contains

  subroutine ResistanceCanopyStomataBallBerry(noahmp, IndexShade)

! ------------------------ Code history -----------------------------------
! Original Noah-MP subroutine: STOMATA
! Original code: Guo-Yue Niu and Noah-MP team (Niu et al. 2011)
! Refactered code: C. He, P. Valayamkunnath, & refactor team (He et al. 2023)
! -------------------------------------------------------------------------

    implicit none

    integer          , intent(in   ) :: IndexShade            ! index for sunlit/shaded (0=sunlit;1=shaded)
    type(noahmp_type), intent(inout) :: noahmp

! local variable
    integer                          :: IndIter               ! iteration index
    integer, parameter               :: NumIter = 3           ! number of iterations
    real(kind=kind_noahmp)           :: RadPhotoActAbsTmp     ! temporary absorbed par for leaves [W/m2]
    real(kind=kind_noahmp)           :: ResistanceStomataTmp  ! temporary leaf stomatal resistance [s/m]
    real(kind=kind_noahmp)           :: PhotosynLeafTmp       ! temporary leaf photosynthesis [umol co2/m2/s]
    real(kind=kind_noahmp)           :: NitrogenFoliageFac    ! foliage nitrogen adjustment factor (0 to 1)
    real(kind=kind_noahmp)           :: CarboxylRateMax       ! maximum rate of carbonylation [umol co2/m2/s]
    real(kind=kind_noahmp)           :: MPE                   ! prevents overflow for division by zero
    real(kind=kind_noahmp)           :: RLB                   ! boundary layer resistance [s m2 / umol]
    real(kind=kind_noahmp)           :: TC                    ! foliage temperature [deg C]
    real(kind=kind_noahmp)           :: CS                    ! co2 concentration at leaf surface [Pa]
    real(kind=kind_noahmp)           :: KC                    ! co2 Michaelis-Menten constant [Pa]
    real(kind=kind_noahmp)           :: KO                    ! o2 Michaelis-Menten constant [Pa]
    real(kind=kind_noahmp)           :: A,B,C,Q               ! intermediate calculations for RS
    real(kind=kind_noahmp)           :: R1,R2                 ! roots for RS
    real(kind=kind_noahmp)           :: PPF                   ! absorb photosynthetic photon flux [umol photons/m2/s]
    real(kind=kind_noahmp)           :: WC                    ! Rubisco limited photosynthesis [umol co2/m2/s]
    real(kind=kind_noahmp)           :: WJ                    ! light limited photosynthesis [umol co2/m2/s]
    real(kind=kind_noahmp)           :: WE                    ! export limited photosynthesis [umol co2/m2/s]
    real(kind=kind_noahmp)           :: CP                    ! co2 compensation point [Pa]
    real(kind=kind_noahmp)           :: CI                    ! internal co2 [Pa]
    real(kind=kind_noahmp)           :: AWC                   ! intermediate calculation for wc
    real(kind=kind_noahmp)           :: J                     ! electron transport [umol co2/m2/s]
    real(kind=kind_noahmp)           :: CEA                   ! constrain ea or else model blows up
    real(kind=kind_noahmp)           :: CF                    ! [s m2/umol] -> [s/m]
    real(kind=kind_noahmp)           :: T                     ! temporary var
! local statement functions
    real(kind=kind_noahmp)           :: F1                    ! generic temperature response (statement function)
    real(kind=kind_noahmp)           :: F2                    ! generic temperature inhibition (statement function)
    real(kind=kind_noahmp)           :: AB                    ! used in statement functions
    real(kind=kind_noahmp)           :: BC                    ! used in statement functions
    F1(AB, BC) = AB**( (BC - 25.0) / 10.0 )
    F2(AB)     = 1.0 + exp( (-2.2e05 + 710.0 * (AB + 273.16)) / (8.314 * (AB + 273.16)) )

! --------------------------------------------------------------------
    associate(                                                                        &
              PressureAirRefHeight    => noahmp%forcing%PressureAirRefHeight         ,& ! in,  air pressure [Pa] at reference height
              TemperatureAirRefHeight => noahmp%forcing%TemperatureAirRefHeight      ,& ! in,  air temperature [K] at reference height
              SoilTranspFacAcc        => noahmp%water%state%SoilTranspFacAcc         ,& ! in,  accumulated soil water transpiration factor (0 to 1)
              IndexGrowSeason         => noahmp%biochem%state%IndexGrowSeason        ,& ! in,  growing season index (0=off, 1=on)
              NitrogenConcFoliage     => noahmp%biochem%state%NitrogenConcFoliage    ,& ! in,  foliage nitrogen concentration [%]
              NitrogenConcFoliageMax  => noahmp%biochem%param%NitrogenConcFoliageMax ,& ! in,  foliage nitrogen concentration when f(n)=1 [%]
              QuantumEfficiency25C    => noahmp%biochem%param%QuantumEfficiency25C   ,& ! in,  quantum efficiency at 25c [umol co2 / umol photon]
              CarboxylRateMax25C      => noahmp%biochem%param%CarboxylRateMax25C     ,& ! in,  maximum rate of carboxylation at 25c [umol co2/m**2/s]
              CarboxylRateMaxQ10      => noahmp%biochem%param%CarboxylRateMaxQ10     ,& ! in,  change in maximum rate of carboxylation for each 10C temp change
              PhotosynPathC3          => noahmp%biochem%param%PhotosynPathC3         ,& ! in,  C3 photosynthetic pathway indicator: 0. = c4, 1. = c3
              SlopeConductToPhotosyn  => noahmp%biochem%param%SlopeConductToPhotosyn ,& ! in,  slope of conductance-to-photosynthesis relationship
              Co2MmConst25C           => noahmp%energy%param%Co2MmConst25C           ,& ! in,  co2 michaelis-menten constant at 25c [Pa]
              O2MmConst25C            => noahmp%energy%param%O2MmConst25C            ,& ! in,  o2 michaelis-menten constant at 25c [Pa]
              Co2MmConstQ10           => noahmp%energy%param%Co2MmConstQ10           ,& ! in,  q10 for Co2MmConst25C
              O2MmConstQ10            => noahmp%energy%param%O2MmConstQ10            ,& ! in,  q10 for ko25
              ConductanceLeafMin      => noahmp%energy%param%ConductanceLeafMin      ,& ! in,  minimum leaf conductance [umol/m**2/s]
              TemperatureCanopy       => noahmp%energy%state%TemperatureCanopy       ,& ! in,  vegetation temperature [K]
              VapPresSatCanopy        => noahmp%energy%state%VapPresSatCanopy        ,& ! in,  canopy saturation vapor pressure at TV [Pa]
              PressureVaporCanAir     => noahmp%energy%state%PressureVaporCanAir     ,& ! in,  canopy air vapor pressure [Pa]
              PressureAtmosO2         => noahmp%energy%state%PressureAtmosO2         ,& ! in,  atmospheric o2 pressure [Pa]
              PressureAtmosCO2        => noahmp%energy%state%PressureAtmosCO2        ,& ! in,  atmospheric co2 pressure [Pa]
              ResistanceLeafBoundary  => noahmp%energy%state%ResistanceLeafBoundary  ,& ! in,  leaf boundary layer resistance [s/m]
              RadPhotoActAbsSunlit    => noahmp%energy%flux%RadPhotoActAbsSunlit     ,& ! in,  average absorbed par for sunlit leaves [W/m2]
              RadPhotoActAbsShade     => noahmp%energy%flux%RadPhotoActAbsShade      ,& ! in,  average absorbed par for shaded leaves [W/m2]
              ResistanceStomataSunlit => noahmp%energy%state%ResistanceStomataSunlit ,& ! out, sunlit leaf stomatal resistance [s/m]
              ResistanceStomataShade  => noahmp%energy%state%ResistanceStomataShade  ,& ! out, shaded leaf stomatal resistance [s/m]
              PhotosynLeafSunlit      => noahmp%biochem%flux%PhotosynLeafSunlit      ,& ! out, sunlit leaf photosynthesis [umol co2/m2/s]
              PhotosynLeafShade       => noahmp%biochem%flux%PhotosynLeafShade        & ! out, shaded leaf photosynthesis [umol co2/m2/s]
             )
! ----------------------------------------------------------------------

    ! initialization
    MPE = 1.0e-6

    ! initialize ResistanceStomata=maximum value and photosynthesis=0 because will only do calculations
    ! for RadPhotoActAbs  > 0, in which case ResistanceStomata <= maximum value and photosynthesis >= 0
    CF = PressureAirRefHeight / (8.314 * TemperatureAirRefHeight) * 1.0e06  ! unit conversion factor
    ResistanceStomataTmp = 1.0 / ConductanceLeafMin * CF
    PhotosynLeafTmp      = 0.0
    if ( IndexShade == 0 ) RadPhotoActAbsTmp = RadPhotoActAbsSunlit  ! Sunlit case
    if ( IndexShade == 1 ) RadPhotoActAbsTmp = RadPhotoActAbsShade   ! Shaded case

    ! only compute when there is radiation absorption
    if ( RadPhotoActAbsTmp > 0.0 ) then

       NitrogenFoliageFac = min(NitrogenConcFoliage/max(MPE, NitrogenConcFoliageMax), 1.0)
       TC                 = TemperatureCanopy - ConstFreezePoint
       PPF                = 4.6 * RadPhotoActAbsTmp
       J                  = PPF * QuantumEfficiency25C
       KC                 = Co2MmConst25C * F1(Co2MmConstQ10, TC)
       KO                 = O2MmConst25C * F1(O2MmConstQ10, TC)
       AWC                = KC * ( 1.0 + PressureAtmosO2 / KO )
       CP                 = 0.5 * KC / KO * PressureAtmosO2 * 0.21
       CarboxylRateMax    = CarboxylRateMax25C / F2(TC) * NitrogenFoliageFac * &
                            SoilTranspFacAcc * F1(CarboxylRateMaxQ10, TC)
       ! first guess ci
       CI  = 0.7 * PressureAtmosCO2 * PhotosynPathC3 + 0.4 * PressureAtmosCO2 * (1.0 - PhotosynPathC3)
       ! ResistanceLeafBoundary: s/m -> s m**2 / umol
       RLB = ResistanceLeafBoundary / CF
       ! constrain PressureVaporCanAir
       CEA = max(0.25*VapPresSatCanopy*PhotosynPathC3 + 0.40*VapPresSatCanopy*(1.0-PhotosynPathC3), &
                 min(PressureVaporCanAir,VapPresSatCanopy))

       ! ci iteration
       do IndIter = 1, NumIter
          WJ = max(CI-CP, 0.0) * J / (CI + 2.0*CP) * PhotosynPathC3 + J * (1.0 - PhotosynPathC3)
          WC = max(CI-CP, 0.0) * CarboxylRateMax / (CI + AWC) * PhotosynPathC3 + &
               CarboxylRateMax * (1.0 - PhotosynPathC3)
          WE = 0.5 * CarboxylRateMax * PhotosynPathC3 + &
               4000.0 * CarboxylRateMax * CI / PressureAirRefHeight * (1.0 - PhotosynPathC3)
          PhotosynLeafTmp = min(WJ, WC, WE) * IndexGrowSeason
          CS = max(PressureAtmosCO2-1.37*RLB*PressureAirRefHeight*PhotosynLeafTmp, MPE)
          A  = SlopeConductToPhotosyn * PhotosynLeafTmp * PressureAirRefHeight * CEA / &
               (CS * VapPresSatCanopy) + ConductanceLeafMin
          B  = (SlopeConductToPhotosyn * PhotosynLeafTmp * PressureAirRefHeight / CS + ConductanceLeafMin) * &
               RLB - 1.0
          C  = -RLB
          if ( B >= 0.0 ) then
             Q = -0.5 * (B + sqrt(B*B-4.0*A*C))
          else
             Q = -0.5 * (B - sqrt(B*B-4.0*A*C))
          endif
          R1   = Q / A
          R2   = C / Q
          ResistanceStomataTmp = max(R1, R2)
          CI   = max(CS-PhotosynLeafTmp*PressureAirRefHeight*1.65*ResistanceStomataTmp, 0.0)
       enddo

       ! ResistanceStomata:  s m**2 / umol -> s/m
       ResistanceStomataTmp = ResistanceStomataTmp * CF

    endif ! RadPhotoActAbsTmp > 0.0

    ! assign updated values
    ! Sunlit case
    if ( IndexShade == 0 ) then
       ResistanceStomataSunlit = ResistanceStomataTmp
       PhotosynLeafSunlit      = PhotosynLeafTmp
    endif
    ! Shaded case
    if ( IndexShade == 1 ) then
       ResistanceStomataShade  = ResistanceStomataTmp
       PhotosynLeafShade       = PhotosynLeafTmp
    endif

    end associate

  end subroutine ResistanceCanopyStomataBallBerry

end module ResistanceCanopyStomataBallBerryMod
